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VORTEX  NUCLEATION  IN  A  DISSIPATIVE  VARIANT  OF  THE 
NONLINEAR  SCHRODINGER  EQUATION  UNDER  ROTATION 

R.  CARRETERO-GONZALEZt,  RG.  KEVREKIDIS*,  AND  T.  KOLOKOLNIKOV§ 

Abstract.  In  the  present  work,  we  motivate  and  explore  the  dynamics  of  a  dissipative  variant 
of  the  nonlinear  Schrodinger  equation  under  the  impact  of  external  rotation.  As  in  the  well  estab¬ 
lished  Hamiltonian  case,  the  rotation  gives  rise  to  the  formation  of  vortices.  We  show,  however,  that 
the  most  unstable  mode  leading  to  this  instability  scales  with  an  appropriate  power  of  the  chemical 
potential  fi  of  the  system,  increasing  proportionally  to  n2'3.  The  precise  form  of  the  relevant  for¬ 
mula,  obtained  through  our  asymptotic  analysis,  provides  the  most  unstable  mode  as  a  function  of 
the  atomic  density  and  the  trap  strength.  We  show  how  these  unstable  modes  typically  nucleate  a 
large  number  of  vortices  in  the  periphery  of  the  atomic  cloud.  However,  through  a  pattern  selection 
mechanism,  prompted  by  symmetry-breaking,  only  few  isolated  vortices  are  pulled  in  sequentially 
from  the  periphery  towards  the  bulk  of  the  cloud  resulting  in  highly  symmetric  stable  vortex  config¬ 
urations  with  far  fewer  vortices  than  the  original  unstable  mode.  These  results  may  be  of  relevance 
to  the  experimentally  tractable  realm  of  finite  temperature  atomic  condensates. 

Key  words.  Vortex  nucleation,  nonlinear  Schrodinger  equation,  Gross-Pitaevskii  equation, 
Bose-Einstein  condensates. 
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1.  Introduction.  Vortices  are  persistent  circulating  flow  patterns  that  occur  in 
many  diverse  scientific  and  mathematical  contexts  [1],  ranging  from  hydrodynamics, 
superfluids,  and  nonlinear  optics  [2,  3]  to  specific  case  examples  in  sunspots,  dust 
devils  [4],  and  plant  propulsion  [5].  The  realm  of  atomic  Bose-Einstein  condensates 
(BECs)  [6,  7,  8]  has  produced  a  novel  and  pristine  setting  where  numerous  features 
of  the  exciting  nonlinear  dynamics  of  single-  and  multi-charge  vortices,  as  well  as  of 
vortex  crystals  and  vortex  lattices,  can  be  not  only  theoretically  studied,  but  also 
experimentally  observed. 

The  first  experimental  observation  of  vortices  in  atomic  BECs  [9]  by  means  of 
a  phase-imprinting  method  between  two  hyperfine  spin  states  of  a  87Rb  BEC  [10] 
paved  the  way  for  a  systematic  investigation  of  their  dynamical  properties.  Stirring 
the  BECs  [11]  above  a  certain  critical  angular  speed  [12,  13,  14]  led  to  the  production 
of  few  vortices  [14]  and  even  of  robust  vortex  lattices  [15,  16].  Other  vortex-generation 
techniques  were  also  used  in  experiments,  including  the  breakup  of  the  BEC  super¬ 
fluidity  by  dragging  obstacles  through  the  condensate  [17],  as  well  as  nonlinear  in¬ 
terference  between  condensate  fragments  [18].  In  addition,  apart  from  unit-charged 
vortices,  higher-charged  vortex  structures  were  produced  [19]  and  their  dynamical 
(in)stability  was  examined.  To  these  earlier  experimental  developments,  one  can  add 
in  recent  years:  the  formation  of  vortices  through  a  quench  of  a  gas  of  atoms  from 
well  above  to  well-below  the  BEC  transition  via  the  so-called  Kibble- Zurek  mecha¬ 
nism  [20];  the  dynamical  visualization  of  such  “nucleated”  vortices  [21]  and  even  of 
vortex  pairs,  i.e.,  dipoles  consisting  of  two  oppositely  charged  vortices;  the  nucleation 
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of  dipoles  via  the  dragging  of  a  laser  beam  through  the  BEC  [22];  the  systematic 
experimental  exploration  of  dipole  dynamics  [23] ;  the  generation  (via  instabilities)  of 
3- vortex  configurations  of  same  or  opposite  signs  in  Ref.  [24]  and  the  “dialing  in”  of 
arbitrary  numbers  of  (few)  same-charge  vortices  and  the  visualization  of  their  intrigu¬ 
ing,  potentially  symmetry-breaking  dynamics  [25].  Naturally,  the  above  developments 
suggest  that  the  study  of  vortices  and  of  their  nucleation  in  BECs  is  a  theme  of  broad 
and  intense  ongoing  interest. 

On  the  other  hand,  another  topic  receiving  increasing  attention  has  concerned 
the  role  of  finite-temperature  induced  “damping”  of  the  BEC  [26].  A  wide  range 
of  recent  examples  has  indicated  that  this  leads  to  anti-damping  motion  of  the  co¬ 
herent  structures  such  as  solitary  waves  (dark  solitons)  and  vortices.  Early  soliton 
experiments  of  about  15  years  ago  observed  the  motion  of  a  dark  soliton  towards 
the  edge  of  the  trap  [27,  28,  29].  It  is  interesting  to  note,  however,  that  this  type 
of  anti-damping  effect  has  been  observed  in  a  far  more  pronounced  way  in  recent 
experiments  of  dark  soliton  oscillations  in  a  unitary  Fermi  gas  [30] .  A  number  of  the¬ 
oretical  studies  have  provided  relevant  explanation  for  this  phenomenology  in  atomic 
BECs  [31,  32,  33,  34,  35,  36,  37,  38,  39].  In  particular,  it  has  been  identified  in  these 
works  that  the  dark  soliton  follows  an  anti-damped  harmonic  oscillator  behavior,  lead¬ 
ing  to  trajectories  of  growing  amplitude  around  the  center  of  the  trap,  until  expelled 
from  the  BEC.  This  motion  has  been  observed  in  the  context  of  the  so-called  dissipa¬ 
tive  Gross-Pitaevskii  equation  model  (DGPE).  The  DGPE  was  originally  introduced 
phenomenologically  by  Pitaevskii  [40]  as  a  way  to  use  a  damping  term  to  account  for 
the  role  of  finite  temperature  induced  fluctuations  in  the  BEC  dynamics;  see,  e.g., 
Refs.  [41,  42,  43,  44]  for  discussions  and  microscopic  interpretations  of  such  a  term. 
Comparisons  [33]  of  its  results  with  more  elaborate  models  such  as  the  (averaged 
quantities  of)  the  stochastic  Gross-Pitaevskii  equation  (SGPE)  [45,  46]  offered  good 
reason  for  exploring  the  simpler  DGPE  model,  as  regards  coherent  structure  (such  as 
soliton)  dynamics. 

More  importantly  for  our  theme  of  vortex  dynamics,  an  increasing  volume  of 
literature  has  been  exploring  the  role  of  thermal  effects  [47,  48,  49,  50,  51].  Here, 
too,  and  in  accordance  with  theoretical  predictions  [52]  (see  also  Ref.  [53]  and  for  a 
recent  discussion  [54]),  an  outward,  in  this  case  spiraling,  trajectory  is  found  for  the 
single  vortex  motion  which  leads  to  its  expulsion  from  the  trap.  While  important 
aspects  of  the  vortex  dynamics  in  the  presence  of  the  thermal  component  such  as 
the  single  vortex  motion  [52,  53]  and  even  the  vortex-pair  interaction  [54]  have  been 
explored  theoretically  (and  numerically),  to  the  best  of  our  knowledge,  the  predictions 
of  the  DGPE  model  in  the  context  of  vortex  nucleation  under  rotation  have  not  been 
previously  examined. 

The  principal  scope  of  our  study  will,  thus,  be  to  provide  some  insight  on  the 
instability  and  dynamics  that  leads  to  the  emergence  of  vortices  in  the  presence  of 
“thermal  dissipation”,  i.e.,  in  the  DGPE  framework.  In  fact,  to  facilitate  the  analy¬ 
sis,  we  will  go  one  step  further  in  simplifying  the  problem  and  will  also  explore  the 
“imaginary  time”  analogue  of  the  GPE.  This  choice  will  be  suitably  explained  and 
motivated  in  the  next  section.  Subsequently,  in  Section  3,  we  will  provide  the  analysis 
that  identifies  the  most  unstable  eigenmode  and  its  scaling  with  the  system  parame¬ 
ters  (most  notably,  the  chemical  potential  /i).  Finally,  in  Section  4,  we  will  summarize 
our  findings  and  present  a  number  of  directions  for  future  study. 

2.  Model  Setup:  the  NLS  Equation  Under  Rotation  and  its  Dissipative 
Variant. 
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2.1.  Dissipationless  case.  The  standard  GPE  model  valid  at  T  =  0  for  describ¬ 
ing  the  the  quasi-2D  condensate  wavefunction  u(x,y,t)  in  the  presence  of  rotation  is: 

iut  =  +  tp?Uapr2u  -  (M+  \u\2u  +  i£lvotue,  (2.1) 

where  (•)*  =  d(-)/dt  and  (■)#  =  d(-)/dQ  and  (r,  0)  are  the  polar  coordinates.  Here 
the  potential  is  assumed  as  representing  a  parabolic  (typically  induced  magnetically) 
trap  of  strength  fltrap,  while  an  external  rotation  of  strength  firot  is  assumed  to  be 
imposed.  We  have  also  explicitly  included  the  chemical  potential  (i  in  the  model 
(although  it  can  be  factored  out  by  a  gauge  transformation),  as  it  will  be  relevant 
in  the  DGPE  variant  of  the  system.  Notice  that  here  we  use  the  dimensionless  form 
of  the  pancake-shaped,  2D  BEC  model  that  has  been  well  established  in  a  variety  of 
archival  references  in  the  field  [7,  8,  55]. 

Before  we  move  to  the  DGPE  variant  of  the  model,  we  should  note  a  remarkable 
implication  of  Eq.  (2.1).  In  particular,  when  analyzing  the  spectrum  of  a  particular 
state  uo,  by  performing  Bogolyubov-de  Gennes  (BdG)  analysis  to  explore  its  stability, 
the  equation  obtained  for  u  =  uq(x,  y )  +  ev(x,  y ,  t)  is  of  the  form: 

ivt  =  +  t^apr2v  -  nv  +  2\u0\2v  +  ufa*  +  iQm tvg,  (2.2) 

where  (•)*  denotes  complex  conjugation.  Now,  decomposing  the  perturbation  as 
v(x,y,t)  =  a{r)eirn6 elujt  +  b{r)e~irn6 e~lujt ,  it  is  straightforward  to  see  that  for  a  ra¬ 
dial  state  (such  as  the  ground  state  of  the  system),  the  sole  influence  of  the  rotation 
frequency  Drot  is  to  shift  the  frequencies  u)  uj  ±  mQro t . 

This  is  illustrated,  e.g.,  in  Fig.  2.1  by  direct  numerical  computations  involving 
the  BdG  linearization  around  the  ground  state  for  the  case  of  fitrap  =  0.2  for  2 
different  values  of  fi  =  5  (left)  and  y  =  10  (right).  The  lowest,  well-known  modes  of 
the  condensate  dynamics  namely  the  dipolar,  quadrupolar,  and  hexapolar  modes  at, 
respectively,  uj  =  Dtrap5  w  =  V^fltrap,  and  uo  =  V^fltrap  (all  with  double  degeneracy), 
are  shown  by  thick  green,  orange  and  yellow  lines  respectively,  showcasing  the  validity 
of  the  above  eigenfrequency  shift  statement.  However,  there  are  numerous  additional 
intriguing  features  to  observe  in  the  figure.  For  one  thing,  we  note  that  since  the 
ground  state  is  stable  and  its  imaginary  eigenvalues  (real  eigenfrequencies)  shift  along 
the  imaginary  axis  of  the  complex  spectral  plane  (Re(A),Im(A)),  the  ground  state 
will  never  become  dynamically  unstable.  Instead,  what  happens  is  that  it  becomes 
energetically  unstable  acquiring  what  is  known  as  negative  energy  modes  [56]  or  in  the 
mathematical  literature  as  negative  Krein  signature  modes  [57].  These  modes  indicate 
that  while  the  solution  may  not  be  dynamically  unstable,  it  is  no  longer  the  ground 
state  of  the  system.  Moreover,  if  a  pathway,  such  as  the  presence  of  dissipation, 
becomes  available  for  relaxing  to  the  ground  state  of  the  system  then  it  would  do 
so  [58]. 

Some  additional  observations  are  also  in  order.  In  particular,  it  is  worthwhile 
to  note  that  among  the  modes  crossing  zero  to  become  negative  energy  or  signature 
ones,  it  is  neither  the  m  =  1,  nor  the  m  —  2  ones  that  do  this  the  first.  Instead, 
modes  associated  with  higher  m  (but  which  start  at  larger  cv)  values  move  faster  and 
cross  0,  spearheading  the  energetic  instability  of  the  present  state.  For  instance,  as 
is  depicted  in  Fig.  2.1,  the  modes  with  m  =  mc  =  15  and  m  =  mc  =  23  are  the  first 
modes  to  cross  the  energetic  stability  threshold  for,  respectively,  y  =  5  and  y  =  10. 
This  instability  occurs  at  firot  =  firot,c  ~  0.349  fltrap  fc>r  /x  =  10  and  at  firot  = 
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flrot/  ^trap 


fhot/  ^trap 


Fig.  2.1.  (color  online)  The  spectrum  of  imaginary  eigenvalues  (normalized  to  the  trap  fre¬ 
quency)  of  the  Hamiltonian  linearization  problem  of  Eq.  (2.2)  in  the  presence  of  rotation.  The  left 
panel  is  for  /a  =  5,  while  the  right  one  is  for  /a  =  10;  £ltrap  =0.2  was  chosen.  The  thick  green, 
orange  and  yellow  lines  correspond  to  the  lowest  three  modes  of  m  —  1,2,3,  giving  the  theoretical 
prediction  of  how  they  depend  as  a  function  of  frequency  according  to  X  =  i(iv(m  =  0)  d=  mOrot). 
The  thick  pink  line  corresponds  to  the  m  =  mc  mode  that  first  becomes  unstable  as  the  rotation  is 
increased.  mc  =  15  for  p  =  5  and  mc  =  23  for  p  =  10. 


f}rot)C  ~  0.28f2trap  for  /i  =  5.  Notice  that  the  critical  rotation  threshold  for  ft  =  5  is 
larger  that  the  one  for  ft  =  10,  an  important  feature  to  which  we  will  return  below. 
The  reason  why  the  above  observations  are  especially  interesting  is  the  following.  As 
proved  rigorously  in  Ref.  [57],  the  inclusion  of  dissipation  in  a  Hamiltonian  model 
leads  modes  of  different  energy  (Krein  signature)  to  move  differently,  due  to  their 
distinct  topological  characteristics.  More  specifically,  modes  with  positive  signature 
move  to  the  left  of  the  spectral  plane  becoming  stable/attracting  eigendirections  for 
the  dynamics.  However,  modes  with  negative  Krein  signature  move  in  the  opposite 
direction  of  the  spectral  plane,  namely  to  the  right  hand  plane,  becoming  immediately 
unstable  as  soon  as  the  dissipation  is  turned  on.  This  statement  goes  hand-in-hand 
with  the  opening  of  relaxation  channels  through  which  the  solution  can  now  revert  to 
its  preferred  ground  state  equilibrium,  given  its  energetic  instability. 

2.2.  Dissipative  Gross-Pitaevskii  Equation.  Now,  let  us  project  the  above 
conclusions  to  the  case  of  the  DGPE  which  is  of  the  form  [40]: 

(i  -  pUt  =  -^Au  +  ig*rapr2'U  -  flu  +  \u\2u  +  iOIotue,  (2.3) 

where  7  (>  0)  refers  to  the  temperature  dependent  parameter  that  has  been  discussed 
extensively  in  this  context  [26,  41,  42,  43,  44].  Exploring  a  nearly  realistic  (although 
slightly  higher  than  relevant,  for  illustration  purposes;  see,  e.g.,  the  discussion  of 
Ref.  [54])  value  of  7  =  0.01,  we  obtain  the  results  for  p  =  10  illustrated  in  Figs.  2.2 
and  2.3.  The  former  one  shows  the  spectral  planes  for  4  values  of  the  ratio  ^rot/^trap? 
two  below,  one  (approximately)  at,  and  one  above  the  energetic  instability  threshold  of 


Vortex  Nucleation  in  a  Dissipative  Variant  of  the  NLS  Equation  under  Rotation 


5 


•'  •  0  1 
••  ••  J 

8  .  0°°“^ 

0 

(b)  fir„t  =  0.2  fitrap 

O 

O 

O 

O 

°  O 

O 

O 

O 

O 

0 

_ ' _ 1 _ , _ _ _  O _ _ 

-0.25  -0.2  -0.15  -0.1  -0.05  0  0.05 

Re(A)/S2trap 


0 

© 

© 

O 

O 

* 

© 

(d)Oot  =0.35fitrap  °  °# 

/ 

00 

0  0 

0  0 

°o° 

Q> 

+  O 

O  O 

O 

O  0 

0 

0 

0  0 

0 

0 

0 

c 9 

©  0 
©  © 

©0 

\ 

’2  -0.06  -0.04  -0.02  0  0.02 

Re(A)/Otrap 


Fig.  2.2.  (color  online)  Spectral  planes  increasing  rotation:  (a)  £2rot  =  0,  (b)  Orot/^trap  =  0.2; 
(c)  f^rot/A^trap  =  0.28  and  (d)  Ahot/^trap  =  0.35,  all  for  7  =  0.01  and  fi  =  10.  Eigenvalues  with 
positive  and  negative  Krein  sign  are  depicted,  respectively,  with  blue  (dark)  and  orange  (light)  points 
and  zero-Krein  eigenvalues  are  depicted  with  black  points.  The  successive  panels  clearly  illustrate  the 
instability  due  to  the  collision  of  opposite  Krein  signature  modes  that  starts  at  flrot,c  =  0.28A2trap 
[see  panel  (c)[. 


the  ground  state.  In  the  presence  of  the  7  term  (irrespectively  of  however  small),  when 
all  the  modes  are  positive  energy  ones,  i.e.,  below  the  threshold  of  Orot?c  =  0.28  f}trap5 
all  of  the  eigenvalues  are  on  the  left  half  plane  [see  Figs.  2.2(a)  and  (b)],  hence  the 
configuration  is  dynamically  stable  as  well  (in  the  DGPE  case).  However ,  above  the 
energetic  instability  threshold  for  the  Hamiltonian  problem,  the  existence  of  negative 
Krein  signature  modes  immediately  leads  to  the  bifurcation  of  unstable  eigenmodes  in 
the  right  half  of  the  spectral  plane  [see  Fig.  2.2(d)]  and  the  configuration  is  dynamically 
unstable  for  the  DGPE.  This  instability  is  manifested  as  a  function  of  f^rot/^trap  f°r 
the  DGPE  case  in  Fig.  2.3  showcasing  that  the  nontrivial  real  parts  of  the  relevant 
eigenvalues  emerge  as  the  threshold  is  crossed.  To  complement  the  instability  picture, 
we  depict  in  Fig.  2.4  the  most  unstable  modes  for  rotations  just  above  the  instability 
threshold.  For  /i  =  5  [see  Fig.  2.4(a)]  the  most  unstable  eigenfunction  for  f^rot/^trap  = 
0.3495  (i.e.,  just  above  the  threshold  £2rot,c/^trap  ~  0.349)  is  the  mode  with  m  —  15 
as  expected  from  Fig.  2.1.  Similarly,  for  fi  =  10  [see  Fig.  2.4(b)],  the  most  unstable 
eigenfunction  for  £2rot/^trap  =  0.2802  (i.e.,  just  above  the  threshold  £2rot,c/^trap  ~ 
0.28)  is  the  mode  with  m  =  23  as  expected  from  Fig.  2.1. 

The  above  stability  considerations  allow  us  to  understand  the  bifurcations  (insta¬ 
bilities)  of  steady  states  bearing  no  initial  vorticity  as  the  rotation  of  the  BEC  cloud 
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Fig.  2.3.  (color  online)  The  largest  eigenvalues  leading  to  the  instability  beyond  the  critical 
value  o/firot/fitrap,  in  £/ie  DGPE  case  with  7  =  0.01. 
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Fig.  2.4.  (color  online)  The  most  unstable  eigenfunction  just  past  the  rotation  threshold.  The 
left  and  right  subpanels  corresponds,  respectively ,  to  the  real  and  imaginary  parts  of  the  most  unstable 
eigenfunction,  (a)  For  fi  =  5  and  Urot/^trap  =  0.3495  >  Orot,c/^trap  ~  0.349  the  most  unstable 
eigenfunction  corresponds  to  m  =  15.  (b)  For  p  =  10  and  Hrot/^trap  =  0.2802  >  Hrot,c/^trap  ~ 
0.28  the  most  unstable  eigenfunction  corresponds  to  m  =  23. 


is  increased.  In  particular,  a  number  Nv  of  vortices  will  be  nucleated  at  the  periphery 
of  the  cloud  through  an  unstable  eigenfunction  invariant  under  rotations  by  2tt/Nv. 
See  for  example  the  eigenfunctions  depicted  in  Fig.  2.4.  However,  it  is  crucial  to  note 
that  this  analysis  only  captures  the  initial  stages  of  the  dynamical  evolution  and  the 
eventual  asymptotic  behavior  may  well  be  different.  This  can  be  due  to  symmetry¬ 
breaking  effects  generated  by  infinitesimally  small,  non- symmetric,  perturbations  that 
will  generically  be  present  in  physical  (and  numerical)  setups.  Therefore,  we  now  ex¬ 
plore  the  dynamics  of  the  above  mentioned  unstable  modes  towards  understanding 
what  they  actually  nucleate  as  the  instability  sets  in.  For  this  purpose  we  have  pro¬ 
duced  long-term  simulations  of  the  DGPE  (2.3)  starting  from  the  stationary  state 
bearing  no  vorticity.  Two  typical  evolutions  are  depicted  in  Figs.  2.5  and  2.6  for,  re¬ 
spectively,  fi  =  5  and  ft  =  10.  We  invite  the  interested  reader  to  see  the  full  movies  at 
this  address:  http://nonlinear.sdsu.edu/  carreter/RotatingBEC.html  [Movies#l  and 
#2].  The  simulations  were  chosen  for  rotations  that  are  slightly  above  critical  so  that 
the  steady  state  with  no  vortices  is  (weakly)  unstable.  Let  us  describe  in  detail  the 
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Fig.  2.5.  (color  online)  Evolution  of  an  initial  steady  state  without  vortices  under  rota¬ 
tion.  The  parameters  are  /a  =  5,  7  =  0.01,  Utrap  =  0.2,  and  Tlrot/Tl  =  0.37.  For  each 
of  the  times  indicated,  we  depict  the  density  (top  sub-rows)  and  phase  (bottom  sub-rows).  The 
windows  for  density  and  phase  are,  respectively,  ( x,y )  E  [—16,16]  X  [—16,16]  and  (x,y)  E 
[—29.5,29.5]  X  [—29.5,29.5].  We  invite  the  interested  reader  to  see  the  full  movie  at  this  address: 
http:// nonlinear.sdsu.edu/  carreter/ Rotating BEC.html  [Movie#!]. 


full  evolution  for  the  first  case  with  fi  —  5  for  Qrot  =  0.37  f2trap  >  ^rot,c  ~  0.349  £2trap 
that  is  depicted  in  Fig.  2.5.  As  it  is  clear  from  the  figure,  for  the  chosen  parameter 
values,  the  steady  state  with  no  vortices  is  unstable  towards  a  mode  with  m  =  17 
vortices  initially  growing  at  the  periphery  of  the  cloud  (see  snapshot  at  t  =  7,000). 
This  mode  is  not  apparent  in  the  density  distribution  since  it  is  outside  the  Thomas- 
Fermi  radius  where  the  density  is  too  low  to  be  able  to  be  picked  up.  However,  the 
phase  distribution  clearly  shows  a  series  of  2tt  windings  that  are  nucleated  at  the 
periphery.  It  is  clear  that  the  growth  of  this  unstable  mode  is  prone  to  asymmetries 
since  it  is  generated  numerically  from  the  noise  inherent  in  the  computation  due  to  its 
finite  precision.  This  effect  would  be  similar  in  the  physical  experiments  where  small 
variations  in  the  initial  density  and  the  trapping  break  the  symmetry  of  the  solution. 
This  asymmetry  is  responsible  for  one  of  the  vortices  to  be  closer  to  the  center  of  the 
cloud  than  its  siblings.  This  selection  mechanism  is  responsible  for  one  of  the  vortices 
to  start  spiralling  inwards  (see  snapshots  t  =  9,  000-11,  000).  It  is  interesting  that,  as 
the  chosen  vortex  rotates  close  to  its  siblings,  it  “pushes”  the  other  vortices  outwards 
and  thus  further  contributes  to  this  selection  mechanism.  After  the  chosen  vortex 
relaxes  at  the  center  of  the  trap,  another  unstable  mode  at  the  periphery  grows  and, 
by  the  same  selection  mechanism  explained  above,  spirals  inwards  (see  snapshots  at 
t  =  12,  000-17,  980).  Then,  the  two  central  vortices  arrange  themselves  into  a  steady 
state  configuration  (see  snapshot  at  t  =  17,000)  with  a  small  perturbation  that  re¬ 
mains  at  the  periphery.  However,  in  this  case,  this  state  is  no  longer  unstable  and 
hence  the  dynamics  is  eventually  attracted  to  it  (see  snapshots  at  t  =  17,  000-20,  000). 
In  fact,  the  resulting  state  with  two  corotating  vortices  in  completely  stable  and  thus 
the  configuration  relaxes  towards  it  and  remains  there. 

A  similar  evolution  is  observed  in  the  case  of  ft  =  10  and  £2rot  =  0.3  £2trap  > 
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Fig.  2.6.  (color  online)  Same  as  in  Fig.  2.5  for  firot/^trap  —  0-3  and  /i  =  10.  The  win¬ 
dows  for  density  and  phase  are ,  respectively,  ( x,y )  G  [—22.5,22.5]  X  [—22.5,22.5]  and  ( x,y )  G 
[—35.8,35.8]  X  [—35.8,35.8].  We  invite  the  interested  reader  to  see  the  full  movie  at  this  address: 
http:// nonlinear.sdsu.edu/  carreter/ Rotating BEC.html  [Movie#2]. 


^rot,c  ~  0.28  ^trap  that  is  depicted  in  Fig.  2.6.  In  this  case,  the  state  with  two 
vortices  in  the  bulk  of  the  condensate  is  still  unstable  and  thus  a  third  vortex  needs 
to  be  pulled  from  the  periphery  inwards  to  finally  create  a  corotating  tripole  that  is 
spectrally  stable. 

The  final  fate  of  the  above  selection  mechanism  can  be,  at  least  in  part,  be  at¬ 
tributed  to  the  fact  that  stationary  corotating  vortex  polygons  with  different  numbers 
of  vortices  have  different  stability  properties  for  fixed  parameter  values.  For  instance, 
in  Fig.  2.7  we  depict  the  steady  states,  their  stability  spectra  and  their  most  unstable 
eigenmode  for  the  case  of  ft  =  5  for  zero,  one  and  two  central  vortices.  As  it  is  evident 
from  the  figure,  the  configurations  bearing  zero  vortices  and  one  vortex  are  unstable 
while  the  configuration  with  two  vortices  is  stable.  This  corroborates  the  dynamical 
evolution  depicted  in  Fig.  2.5  where  the  initial  state  with  zero  vortices  destabilizes 
towards  a  transient  state  with  one  vortex  that,  in  turn,  destabilizes  towards  the  final, 
stable  steady  state  with  two  vortices.  A  similar  stability  analysis  for  fa  =  10  (results 
not  shown  here)  indicates  that  indeed  the  steady  states  with  zero,  one  and  two  vortices 
are  all  unstable,  while  the  state  with  three  is  perfectly  stable;  corroborating  what  is 
observed  in  the  dynamical  evolution  depicted  in  Fig.  2.6. 

The  above  results  prompt  the  question  of  stability  for  configurations  bearing  an 
increasing  number  of  vortices.  For  instance,  as  previously  described,  the  configuration 
without  vortices  destabilizes  when  the  rotation  is  increased,  while  configurations  with 
one  or  more  vortices  become  stable.  However,  polygonal  configurations  have  a  limit 
to  the  number  of  vortices  that  they  can  hold  before  becoming  unstable  (see  Ref.  [59] 
and  references  therein).  This  is  depicted  in  Fig.  2.8  where  the  stability  for  polygonal 
vortex  states  with  3,  4  and  5  vortices  for  the  case  ft  =  5  is  examined.  As  it  is  clear 
from  these  results,  the  polygonal  configurations  with  3  and  4  vortices  are  indeed  stable 
while  the  one  with  5  vortices  is  unstable.  Interestingly,  the  instability  responsible  for 
the  breakup  of  the  5-vortex  configuration  is  not  an  angular  mode  as  in  all  the  cases 
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Fig.  2.7.  (color  online)  Stability  of  steady  states  with  different  number  of  vortices  for  n  =  5, 
7  =  0.01;  £2trap  =  0.2;  and  QTOt/Q  =  0.37.  (a)-(c)  Stability  spectra  for  configurations  with  0,  1  and  2 
vortices  respectively,  (d)-(f)  Corresponding  density  (left  subpanels)  and  phases  (right  subpanels)  for 
these  configurations.  (g),(h)  Most  unstable  eigenfunctions  for  configurations  with  0  and  1  vortices, 
respectively.  The  left  and  right  subpanels  corresponds,  respectively,  to  the  real  and  imaginary  parts 
of  the  most  unstable  eigenfunction. 


presented  previously.  In  fact,  it  is  clear  that  all  the  angular  modes  are  stable  as 
it  can  be  seen  from  the  zoomed-in  version  of  the  spectrum  depicted  in  panel  (d)  of 
the  figure  where  only  the  eigenvalues  associated  with  angular  modes  are  depicted. 
Nonetheless,  the  5-vortex  state  is  indeed  unstable  as  is  evident  by  the  small  cluster 
of  unstable  eigenvalues  enclosed  in  the  small  circle  in  panel  (c)  of  the  figure.  The 
eigenmodes  associated  with  these  unstable  eigenvalues  are  depicted  in  panels  (h)- 
(j).  These  symmetry  breaking  modes  bring  some  of  the  vortices  closer  and  others 
further  apart  from  each  other.  It  is  precisely  this  mechanism  that  is  responsible  for 
the  destabilization  of  the  5-vortex  polygonal  state  as  it  is  depicted  in  the  snapshots 
of  its  evolution  in  Fig.  2.9.  We  invite  the  interested  reader  to  see  the  full  movie 
at  this  address:  http://nonlinear.sdsu.edu/  carreter/RotatingBEC.html  [Movie#3]. 
Here,  the  initial  steady  state  bearing  a  5-vortex  polygonal  state  destabilizes  around 
t  =  8200,  via  a  mode  that  pushes  some  vortices  inward  and  other  outward,  eventually 
resulting  in  a  stable  4-vortex  polygonal  state  after  one  vortex  is  ejected  towards  the 
periphery  of  the  cloud. 

We  have  checked  that  the  above  phenomenology  persists  for  other  values  of  the 
dissipation  coefficient  7.  For  instance,  although  the  stability  thresholds  and  the  order 
(m)  of  the  unstable  modes  for  the  vortex- less  configuration  are  independent  of  7, 
the  growth  rates  for  these  instabilities  are  indeed  dependent  on  dissipation  (see  also 
discussion  below  in  Sec.  3.1).  In  particular,  the  larger  7  is,  the  larger  the  instability 
growth  rates  will  be.  Therefore,  for  larger  values  of  7  the  instability  will  set  in  earlier 
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Fig.  2.8.  (color  online)  Stability  of  steady  states  with  different  number  of  vortices  for  p  =  5, 
7  =  0.01;  12trap  =  0.2;  and  Orot/lhrap  —  0.37.  (a)-(c)  Stability  spectra  for  configurations  with 
3,  4  and  5  vortices  respectively.  Note  the  unstable  eigenvalues  inside  the  circle  for  5  vortices  in 
panel  (c).  (d)  Zoomed-in  version  of  the  stability  spectrum  for  5  vortices  showing  the  eigenvalues 
corresponding  to  angular  modes,  (e)-(g)  Density  (left  subpanels)  and  phases  (right  subpanels)  for 
the  configurations  with  3,  4  and  5  vortices  respectively .  (h)-(j)  Three  most  unstable  eigenfunctions 
for  the  configuration  with  5  vortices.  The  left  and  right  subpanels  corresponds,  respectively ,  to  the 
real  and  imaginary  parts  of  the  most  unstable  eigenfunction. 
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Fig.  2.9.  (color  online)  Evolution  of  an  unstable  5-vortex  configuration.  The  initial  state 
corresponds  to  the  one  in  Fig.  2.8(g)  (namely  Ihot/lhrap  =  0.37,  p  =  5,  7  =  0.01  and  12trap  = 
0.2,).  The  windows  for  density  and  phase  are,  respectively,  {x,y)  E  [—17.5,17.5]  X  [—17.5,17.5]  and 
(x,y)  E  [—19,19]  X  [—19,19].  We  invite  the  interested  reader  to  see  the  full  movie  at  this  address: 
http:// nonlinear.sdsu.edu/  carreter/ Rotating BEC.html  [Movie#3]. 


and,  more  importantly,  the  pattern  selection  mechanism  for  the  setlling  of  a  cluster  of 
vortices  at  the  center  of  the  cloud  will  be  different.  This  is  a  direct  consequence  from 
the  fact  that  the  spiraling  experienced  by  a  vortex  due  to  dissipation  has  a  faster  radial 
rate  as  7  is  increased  [54] .  A  key  effect  of  the  slow  down  in  the  radial  spiraling  rate  as  7 
is  decreased  is  that  the  pattern  selection  mechanism  has  “more  time”  to  select  votices 
and  thus  a  smaller  number  of  vortices  is  eventually  pulled  in  from  the  periphery. 
This  is  precisely  what  we  observe  in  numerical  simulations  where  smaller  values  of  7 
give  rise  to  final  configurations  with  a  smaller  number  of  vortices  (results  not  shown 
here).  For  instance,  for  (i  =  5  we  find  that  for  values  of  7  of  0.5,  0.2,  0.1,  and  0.01, 
a  vortex-less  configuration  evolves  towards  a  stable  steady  state  configuration  with, 
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Fig.  2.10.  (color  online)  Steady  state  configurations  with  N  +  l  vortices  for  12rot/^trap  =  0.37, 
/a  =  5,  7  =  0.01,  and  Qtrap  =  0.2.  (a)  4  +  1  configuration,  (b)  5+1  configuration,  and  (c)  6+1 

configuration.  All  these  configurations  are  stable  for  the  chosen  parameter  values. 


respectively,  9,  7,  3,  and  2  vortices.  The  same  setup  but  for  /i  =  10  yields,  respectively, 
15,  11,  5,  and  3  vortices.  We  invite  the  interested  reader  to  see  the  full  movies  for  all 
of  these  cases  at  this  address:  http://nonlinear.sdsu.edu/  carreter /Rot at ingBEC.html 
[Movies#4-ll]. 

Finally,  let  us  briefly  touch  upon  the  existence  of  other  relevant  vortex  config¬ 
urations.  The  fact  that  polygonal  configurations  with  large  number  of  vortices  lose 
stability  prompts  the  important  question:  what  are  the  remaining  stable  configura¬ 
tions  of  the  system?  For  instance,  in  the  absence  of  rotation,  it  has  been  shown 
that  polygonal  configurations  become  destabilized  towards  asymmetric  configurations 
in  a  symmetry-breaking  pitchfork  bifurcation  [60,  59]  that  has  been  observed  in  ac¬ 
tual  BEC  experiments  [25].  This  instability  occurs  when  the  polygonal  configuration 
increases  its  radius,  namely,  when  the  angular  momentum  of  the  vortex  cluster  is  in¬ 
creased.  In  a  similar  manner,  as  we  increase  the  rotation  in  the  DGPE  model  (2.3)  or 
as  we  increase  the  number  of  vortices,  polygonal  configurations  lose  stability.  These 
instabilities  can  be  manifested  through  angular  modes  (cf.  the  case  of  one  vortex  in 
Fig.  2.7)  or  through  symmetry  breaking  modes  (cf.  the  case  for  5  vortices  in  Figs.  2.8 
and  2.9).  However,  on  the  other  hand,  if  ones  starts  with  a  polygonal  state  with 
an  extra  vortex  at  the  center,  the  so-called  N+l  vortex  configurations,  new  stable 
states  are  produced  [59]  that  might  even  be  the  ground  states  (i.e.,  the  minimal  en¬ 
ergy  states)  of  the  system  [60].  In  Fig.  2.10  we  depict  three  examples  of  these  N+l 
configurations  for  4+1,  5+1,  and  6+1  vortices.  It  is  important  to  mention  that  these 
three  configuration  are  indeed  stable  for  the  chosen  parameter  values.  In  fact,  as 
the  rotation  increases  and  the  number  of  vortices  increases  as  well,  more  complex 
configurations  relating  to  triangular  (Abrikosov)  vortex  lattices  arise. 

Although  our  principal  emphasis  in  this  work  lies  in  the  identification  of  the 
most  unstable  mode  that  results  in  the  instability  of  the  state  bearing  no  vortices  in 
the  context  of  Eq.  (2.3),  clearly  numerous  additional  issues  emerge  from  the  above 
simulations.  These  concern  the  identification  of  the  minimal  energy  state  and  the 
transition  pathways  resulting  in  different  types  of  local  or  global  attractors.  We  will 
briefly  return  to  these  questions  in  the  context  of  future  challenges  in  Section  4. 

3.  Overdamped  NLS  Model. 

3.1.  Model  and  Numerical  Results.  Now  a  crucial  observation  allows  us  to 
explore  the  principal  question  raised  above  (about  the  dominant  instability  mode) 
in  the  context  of  slightly  different  and  somewhat  simpler  model.  The  topological 
nature  of  the  observed  stability  characteristics  (i.e.,  of  positive  and  negative  energy 
modes)  of  the  stationary  state  without  vortices  renders  them  robust  and  independent 
of  the  precise  value  of  7.  In  particular,  for  different  values  of  7,  the  actual  size  of 
the  growth  rates  will  change  (in  fact,  the  higher  7  is,  the  more  unstable  the  relevant 
modes  are).  However,  as  can  be  also  numerically  checked,  the  ordering  of  the  relevant 
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0.0000  0.1111  0.11357  0.11503  0.18038  0.25000 


0.25000  0.17835  0.14069  0.10147  0.07180  0.01718 


Fig.  3.1.  (color  online)  Top:  snapshots  of  the  simulation  of  (3.1)  with  slowly  varying  f2rot  (as 
indicated)  and  with  other  parameters  given  by  12trap  =  0.3538  and  /a  =  16.1  (cf.  Ref.  [25]).  Top 
row:  Hrot  is  slowly  increased  from  0  to  0.25  according  to  the  formula  (3.2).  Middle  row:  Hrot  is 
slowly  decreased  from  0.25  to  0  according  to  the  formula  (3.3).  Bottom  row:  number  of  vortices  as 
a  function  of  Hrot  • 


eigenmodes  will  not  be  modified  by  the  precise  value  of  7.  That  is  to  say,  the  same 
mode  will  become  unstable  at  the  same  critical  point  (but  with  a  different  slope  of 
Re(A)  vs.  flrot/^tmp  in  Fig.  2.3)  for  a  different  value  of  7.  Given  this  feature,  we 
will  hereafter  choose  to  explore  the  “over damped”  limit  of  large  7,  where  in  fact  the 
Hamiltonian  term  in  the  left  hand  side  of  Eq.  (2.3)  is  neglected  in  comparison  to 
the  7-dependent  dissipative  one  (and  subsequently  a  time  rescaling  to  absorb  7  is 
performed).  Hence,  we  will  work  below  with  the  “imaginary  time”  variant  of  the 
equation 

ut  =  u  -  ^CL2Uavx2u  +  HU  -  \u\2u  -  ittrotUe,  (3.1) 

which,  based  on  the  above  arguments,  should  be  sufficient  to  provide  us  with  a  pre¬ 
diction  for  the  above  instability  when  focusing  on  the  vortex-less  stationary  state.  In 
fact,  this  very  statement  will  be  double  checked  a  posteriori  in  the  next  section. 

For  motivational  purposes,  we  start  our  presentation  with  two  direct  numerical 
experiments  of  Eq.  (3.1). 

Experiment  1:  Slow  Increase  in  Rotation.  Let  us  allow  for  the  rotation  Qrot 
to  slowly  increase  from  0  to  0.25.  More  specifically,  we  start  with  Qrot  =  0,  evolve 
the  system  until  t  =  2000  (so  that  it  effectively  reaches  its  ground  state  under  the 
imaginary  time  integration),  then  start  increasing  Hrot  linearly  from  0  to  0.25  as  t 
varies  from  2000  to  4000,  and  then  finally  we  keep  f2rot  =  0.25  until  t  =  10000.  In 
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full,  this  experimental  protocol  can  be  summarized  as: 

rot  =  0 . 25*min(max(  (t-2000)  /2000 ,0)  ,  1). 


(3.2) 


We  used  FlexPDE  to  simulate  Eq.  (3.1)  with  a  uniform  mesh  with  around  13000 
triangles  and  adaptive  time  stepping. 

Experiment  2:  Slow  Decrease  in  Rotation  and  Hysteresis.  Let  us  now  slowly 
decrease  Hrot  from  0.25  to  0.  In  this  complementary  numerical  experiment,  we  start 
with  Hrot  =  0.25,  and  dynamically  evolve  the  system  until  t  =  2000  to  allow  it  to 
relax  to  its  preferred  vortex-lattice  ground  state  profile.  Then,  we  start  decreasing 
Hrot  linearly  to  0  as  t  varies  from  2000  to  4000,  then  keep  Hrot  =  0  until  t  =  10000. 
Mathematically  again,  this  procedure  can  be  summarized  as: 


(3.3) 


Hrot  =  0. 25* (l-min(max( (t-2000) /2000, 0) ,  Ill- 


Figure  3.1  presents  a  series  of  snapshots  for  each  of  the  described  experiments.  There 
are  a  number  of  interesting  observations  to  make  here,  as  well  as  connections  to  pro¬ 
vide  with  the  discussion  in  the  previous  section.  As  flrot  is  ramped  up,  up  to  a  critical 
rotational  frequency  no  vortices  arise.  However,  when  they  do  arise  (and  despite  the 
weak  ramping),  a  considerable  number  of  vortices  seems  to  emerge  asymmetrically  (at 
first)  yet  nearly  simultaneously.  Gradually,  as  the  rotation  frequency  increases,  addi¬ 
tional  vorticity  is  “elicited”  from  the  boundary,  eventually  leading  the  configuration  to 
self-organize  into  a  triangular  vortex  lattice  as  the  final  rotation  frequency  is  reached. 
On  the  other  hand,  while  the  frequency  is  decreased  towards  Hrot  =  0,  we  can  see 
that  the  process  is  clearly  hysteretic,  as  for  similar  values  of  the  rotation  frequency 
as  in  the  top  panel,  a  considerably  larger  number  of  vortices  appears  to  survive.  This 
feature  is  most  dramatic  near  the  Hrot  =  0  (e.g.  in  the  next  to  last  panel  of  the 
bottom  row)  where  numerous  vortices  appear  to  survive  in  this  metastable  dynamics, 
although  clearly  this  is  far  from  the  ground  state  for  that  rotation  frequency.  Al¬ 
though  there  are  numerous  features  that  one  may  wish  to  explore  on  the  basis  of  this 
dynamical  simulation  and  the  numerical  results  presented  in  the  previous  section,  the 
one  that  we  will  focus  on  below  (following  up  on  the  discussion  of  the  earliest  part  of 
the  previous  section)  concerns  the  instability  dynamics  of  the  dissipative  variant  of  the 
GPE  for  the  vortex-less  steady  state  configuration.  In  particular,  our  expectation  on 
the  basis  of  the  above  direct  simulation,  as  well  as  from  the  stability  results  presented 
is  that  a  large  m  mode  is  the  one  (predominantly)  responsible  for  the  destabilization 
of  the  vortex- less  state,  as  is  indeed  observed  in  Fig.  3.1.  We  now  proceed  to  analyze 
this  trait  mathematically  in  more  detail  at  the  level  of  Eq.  (3.1). 

3.2.  Asymptotic  Analysis.  To  start  our  analysis,  we  rescale  Eq.  (3.1)  as  fol¬ 
lows: 


After  dropping  the  hats,  we  obtain 


(3.4) 


where 


(3.5) 
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This  is  equivalent  to  the  following  real  system  of  PDEs  for  the  real  and  imaginary 
parts  of  u  =  v  +  iw: 


{Vt  =  e2Av  +  ^1  —  \x\ v  —  v3  —  W2V  +  Qwq, 
wt  =  e2Aw  +  ^1  —  \x\2^j  w  —  w3  —  n2re  —  Qv$. 


Let  u  =  770  be  the  radially-symmetric  vortex- less  state  which  satisfies: 

0  =  £2  ^T/Orr  +  r ^  +  (l  ~  r2)  7?0  -  T]q- 


(3.6) 


(3.7) 


We  linearize  Eq.  (3.6)  around  770,  so  that 


V  =  VO  +  ext(p(r), 
w  =  0  +  extip(r), 


to  obtain  the  system 


|  A cj)  =  s2A(j)  +  (l  -  r2  -  3?7o)  0  + 

\  \i/j  =  s2 Aip  +  (l  —  r2  —  ?7o)  V’  — 

We  now  use  the  radial-polar  decomposition  for  the  perturbations  of  the  form: 
</>(r,  6)  =  eirn9 0(r);  0)  =  eirn0 7/1  (r) , 


(3.8) 


to  obtain 


A 4>  =  £2  (j>rr  +  ~  ^T^>)  +  (l  —  T2  -  3 T%)  <j>  +  imQ.'lp, 

\lp  =  £2  [i>rr  +  \lpr  ~  ^V’)  +  (l  -  T2  -  rfc)  Ip  -  imtt(p, 


Changing  variables  to  i/j  =  ix/j,  and  dropping  the  hat  we  then  get  a  purely  real  system 
A <p  =  £2  (<Prr  +  \<j>r  -  pr</>)  +  (1  -  r2  -  3t7q)  4>  -  mflip, 

<  \  2  '  (3-10) 

A  ip  =  £2  yiprr  +  f'i/’r  -  yripj  +  (l  -  r2  -  rj%)  ip  -  mQ<p. 

It  is  known  that  A  <  0  when  Q  =  0  and  A  >  0  for  sufficiently  large  fi,  as  corroborated 
also  by  our  numerical  computations  in  the  previous  section.  We  therefore  seek  the 
instability  threshold  value  for  Q  for  which  A  =  0.  Thus,  setting  A  =  0  in  the  eigenvalue 
problem,  we  obtain  a  modified  eigenvalue  problem  for  mfl  of  the  form: 


I"  rntli/j  =  £2  (<£rr  +  \<pr  -  +  (1  -  r2  -  3^)  cp, 

|  mQcp  =  e2  (iprr  +  \i>r  ~  ^rV’)  +  (1  -  r2  -  rfc)  ip. 


(3.11) 


To  obtain  an  intuitive  understanding  of  the  situation  ahead  of  the  detailed  analysis, 
let  us  solve  this  problem  numerically  and  then  plot  the  graph  of  m  vs.  ft.  This  is 
shown  in  Fig.  3.2(a)  for  e  =  0.04  and  e  =  0.02.  Recall  that  the  large  density/large 
chemical  potential  limit  is  associated  with  e  0,  hence  the  choice  of  suitably  small 
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Fig.  3.2.  (color  online)  (a)  Solution  to  Eq.  (3.11)  showing  Q  as  a  function  of  m.  For  a  given 
mode  m  on  the  horizontal  axis,  the  vertical  axis  shows  the  threshold  value  of  Q  for  which  this  mode 
becomes  unstable.  The  minimum  of  this  graph  is  the  overall  threshold  where  the  instability  first  sets 
in  as  Q  is  increased  from  zero,  (b)  The  profile  ofpo,  its  Thomas- Fermi  asymptotic  approximation 
go  ~  max((l  —  r2),  0)1/2 ,  as  well  as  the  profile  of  the  eigenfunctions  corresponding  to  the  problem 
(3.11).  (c)  Solution  to  the  reduced  system  (3.13)  as  compared  with  the  full  system  (3.11) 


e.  We  find  that  this  graph  has  a  minimum  which  corresponds  to  the  smallest  value 
of  Q  =  Qc  for  which  the  instability  first  manifests  itself.  Critically,  for  our  discussion, 
this  minimum  depends  on  e.  Our  goal  is  to  characterize  this  minimum  analytically,  as 
well  as  to  compute  the  corresponding  wave  number  m  =  mc  which  should  approximate 
the  instability  eigenmode  that  manifests  itself  as  Q  increases  and  first  crosses  past 
Figure  3.2(b)  shows  the  plot  of  the  leading  eigenfunction  for  e  =  0.02  with  m  —  15 
(corresponding  to  the  critical  threshold  =  0.01383).  The  relevant  mode  appears 
to  be  localized  near  r  ~  1.  To  capture  this,  we  first  have  to  resolve  the  corner  layer 
of  the  steady  state  770  accurately  near  r  =  1.  To  that  effect,  we  zoom  in  at  r  =  1  and 
rescale  according  to: 


r  =  1  +  e2/3y,  r?o  =  e1/3U. 
To  leading  order  then  we  get  from  Eq.  (3.7) 

Uyy  =  2yU  +  U3, 


(3.12) 


which  is  a  rescaled  Painleve  II  transcendent.  It  is  well-known  [61]  that  Eq.  (3.12) 
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admits  a  unique  solution  of  the  form 

f  C  Ai  {y/2y)  as  y  oo 
U~  < 


[  \J —‘ly  as  y  — oo 

where  Ai  is  the  Airy  function  and  C  is  a  constant.  Next,  we  use  the  same  change  of 
variables  in  the  eigenvalue  problem  (3.11).  We  obtain  then  the  reduced  problem 


r  m0Q,0ip  =  4>yy  -  ml4>  +  (—2 y  -  3U 2)  <t>, 

{  moao0  =  -ipyy  -  mlip  +  (—2 y  -  U 2)  ip, 
subject  to  the  boundary  conditions  {0,^}  0  as  \y\  — )>  oo,  where 

m  =  £_2/3rao,  =  54/3Ho- 


(3.13) 


(3.14) 


The  problem  (3.13)  is  solved  numerically.  The  solution  is  shown  in  Fig.  3.2(c). 
Superimposed  are  the  solutions  to  the  full  eigenvalue  problem  (3.11)  for  e  =  0.02  and 
5  =  0.04.  It  can  thus  be  clearly  observed  that  the  scaling  (3.14)  is  indeed  correct.  From 
Fig.  3.2,  the  minimum  is  attained  at  around  Ho,c  ~  2.5  (using  the  5  =  0.02  curve). 
We  now  state  this  result  (a  numerically  assisted  proof  for  the  result  in  provided  in 
the  appendix): 

Main  Result.  There  exist  constants  Ho,c  and  mo,c  whose  approximate  values 
are 


120,c  «  2.529,  m0,c  «  1.111 
such  that  the  following  is  true.  Let 

Llc  —  £  ^  ^2o,c5  rnc  —  s  ^  R^o,c* 

Then  the  vortex-less  steady  state  (3.7)  of  Eq.  (3.4)  is  stable  when  Ll  <  Qc  and  be¬ 
comes  unstable  as  Q  crosses  Llc.  The  fastest- growing  unstable  mode  corresponds  to 
the  oscillations  of  the  boundary  with  the  mode  mc. 

Note  that  in  terms  of  the  original  variables  of  Eq.  (3.1),  the  critical  threshold  is 
given  by 

ftrot,c  =  lV2_4/V1/3^tr/ap  »  1.0036A4-1/3fi?5»;  (3.15) 

mc  =  m0,c22/3M2/3Ht_rap3  *  1.7Q/i2/3n-%3.  (3.16) 

In  Experiment  1  of  the  previous  section  depicted  in  Fig.  3.1,  we  had  y  =  16.1 
and  fltmp  =  0.3538;  these  numbers  correspond  to  the  experimentally  relevant  setting 
of  the  work  of  Ref.  [25].  This  yields  firotjC  ~  0.10  and  mc  ~  22.11.  This  is  in 
excellent  agreement  with  the  actual  numerical  simulations.  In  the  top  row  of  the 
figure,  the  instability  becomes  apparent  around  firot  ~  0.11.  The  instability  actually 
sets  in  shortly  prior  to  this,  but  it  takes  some  time  for  it  to  fully  mature.  Once  the 
instability  is  fully  developed  (fourth  snapshot  from  the  left),  it  results  in  24  vortices, 
in  fair  agreement  with  the  predicted  value  of  mc  «  22  (it  should  be  kept  in  mind  that 
some  of  these  vortices  may  be  in  the  periphery  of  the  cloud  and  hence  may  not  be 
discernible  in  the  density  profile  shown). 

It  is  interesting  finally  to  connect  the  results,  e.g.,  with  those  of  Figs.  2. 1-2.3.  In 
particular,  as  indicated  in  these  results  the  left  panel,  e.g.,  of  Fig.  2.1  corresponds  to 


Vortex  Nucleation  in  a  Dissipative  Variant  of  the  NLS  Equation  under  Rotation 


IT 


li  —  5,  while  the  right  one  to  /i  =  10.  According  to  the  above  scaling  in  the  former 
case  firot.c  ~  0.343,  while  in  the  latter  case  firot,c  ~  0.272.  It  can  be  seen  that  these 
critical  points  are  in  excellent  agreement  with  the  crossing  points  (from  positive  to 
negative  energy  modes)  of  the  Hamiltonian  system  in  Fig.  2.1,  which  involves  the  case 
of  7  =  0.  They  are  also  in  excellent  agreement  with  the  stability  thresholds  of  the 
dissipative  system  used  in  Figs.  2.2  and  2.3,  although  the  latter  only  use  7  =  0.01. 
This  comparison  is  given  to  also,  a  posteriori,  justify  the  use  of  the  analysis  in  the 
framework  of  the  overdamped  limit  of  Eq.  (3.1)  in  the  present  context. 

4.  Conclusions  and  Future  Challenges.  In  summary,  in  the  present  work, 
we  have  explored  the  instability  in  the  presence  of  rotation  of  a  dissipative  variant 
of  the  GPE.  We  have  connected  this  instability  to  the  emergence  of  negative  energy 
(or  Krein  signature)  modes  and  the  phenomenon  of  energetic  (but  not  dynamical) 
instability  of  the  radial  vortex-less  profile  in  the  corresponding  Hamiltonian  system 
in  the  absence  of  dissipation.  We  have  also  connected  it  to  the  emergence  of  a  real 
eigenvalue  corresponding  to  suitably  large  m  (i.e.,  azimuthal  order)  for  the  dissipative 
system.  Moreover,  we  have  argued  that  this  m  should  be  independent  of  7,  but  should 
depend  on  the  trapping  frequency  and  on  the  chemical  potential  (i.e.,  the  maximal 
atomic  density)  of  the  system.  We  have  systematically  developed  a  scaling  law  that 
provides  the  critical  rotation  frequency  as  a  function  of  these  parameters.  We  have 
confirmed  the  connection  of  this  scaling  with  (a)  the  imaginary  time  ( “overdamped” ) 
model  used;  (b)  the  original  Hamiltonian  model  and  (c)  the  intermediate  between 
the  two  dissipative  Gross-Pitaevskii  equation  (DGPE)  model.  Using  direct  numerical 
simulations  we  have  corroborated  the  prediction  of  our  asymptotic  analysis  and  ob¬ 
served  that  these  unstable  modes  with  high  mode  number  m  indeed  nucleate  a  large 
number  of  vortices  at  the  periphery  of  the  atomic  cloud.  However,  we  have  found  that 
despite  the  large  number  of  vortices  at  the  periphery,  for  the  small  values  of  7  chiefly 
considered  herein  in  the  DGPE  (yet  somewhat  larger  than  the  ones  that  have  been 
claimed  as  relevant  for  realistic  experimental  settings;  see  for  a  recent  discussion  [54]) 
few  isolated  vortices  are  pulled  in  sequentially  towards  the  center  of  the  cloud.  The 
process  whereby  single  vortices  are  singled  out  from  the  multi-vortex  “necklace”  at 
the  periphery  is  a  pattern  selection  mechanism  based  on  symmetry-breaking.  This 
sequential  recruiting  of  peripheral  vortices  towards  the  cloud  center  saturates  when 
reaching  a  highly  symmetric  configuration  with  a  few  vortices  at  the  center  of  the 
cloud  that  is  stable  for  the  chosen  parameters  (chemical  potential  and  rotation  rate 
normalized  by  trap  strength). 

These  results  open  a  number  of  interesting  directions  for  further  exploration. 
Admittedly,  our  approach  to  the  vortex  nucleation  problem  is  rather  different  from 
that  of  earlier  works;  see  e.g.  [62]  and  the  relevant  discussion  of  Ref.  [8].  Nevertheless, 
it  would  be  particularly  interesting  to  explore  in  current  experimental  settings  (such 
as  e.g.  [25])  whether  a  value  of  7  can  be  “inferred”  (e.g.  from  the  spiraling  motion 
of  a  vortex;  see  e.g.  the  relevant  discussion  of  [54]).  Based  on  such  a  value,  our 
analysis  and  computation  could  provide  diagnostics  both  for  which  eigenmode  will 
cause  the  instability  of  the  vortex-less  state  and  for  which  asymptotic  state  may  be 
experimentally  observed. 

Our  observations  also  raise  a  related  problem.  Clearly,  for  different  values  of 
7  ranging  from  the  under  damped  limit  of  Figs.  2. 5-2. 9  to  the  over  damped  one  of 
Fig.  3.1,  (for  same  trap  strength  and  chemical  potential)  we  have  conclusively  argued 
that  the  same  eigenmode  and  the  same  critical  frequency  are  generically  responsible  for 
the  observed  instability.  Yet,  the  instability  manifestation  has  dramatically  different 
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outcomes  for  the  different  limits.  This  naturally  prompts  the  question:  what  is  the 
favored  asymptotic  state  (depending  on  the  value  of  7)  and  how  do  we  get  there? 
The  first  and  perhaps  simpler  question  (that  has  been  previously  considered;  see 
e.g.  [63]  for  an  early  example)  is  presumably  one  of  energetic  comparisons  of  the 
states  containing  different  numbers  of  vortices  (incorporating  their  angular  momentum 
contributions).  The  second  and,  arguably,  more  difficult  question  is  one  of  transition 
state  pathways  that  enable  the  nucleation  of  different  multi- vortex  configurations. 
The  latter  may  be  particularly  worth  exploring,  especially  since  our  results  indicate 
that  they  will  be  dependent  on  features  such  as  the  thermal  coupling  parameter  7. 

Appendix.  A  Numerically  Assisted  Proof  Of  The  Main  Result. 

Define  the  operators 

r  L^)  :=  cf>yy  +  (-22,  -  3U2)  0, 

1  L2(<l>)  ■■=  4>yy  +  (-22/  -  u2)  0. 

As  the  first  step,  we  show  that  both  L\  and  L2  are  negative  operators.  First,  note 
that  that  both  are  self-adjoint  so  it  suffices  to  show  that  all  eigenvalues  are  real  non¬ 
positive.  To  show  that  L2  is  negative,  simply  note  that  L,2(U)  =  0.  Since  V  is 
positive,  Sturm’s  eigenvalue  oscillation  theorem  then  implies  that  A  =  0  is  the  largest 
eigenvalue  of  L2,  hence  L2  is  negative.  The  negativity  of  L\  follows  from  the  fact  that 
2 y  +  3 U2  >  0  for  all  y\  see  Ref.  [64]. 

The  problem  (3.13)  may  then  be  reformulated  as 

m0O 0  =  ±Vm;  =  (Li  -  ml)  •  (1/2  -  ml)  (f>.  (A.l) 

It  follows  by  the  negativity  of  L\  —  ml  and  L-2  —  ml  that  /j  is  real  and  positive  so  that 
the  curve  Do  =  Do  (mo)  is  well  defined. 

Finally,  we  show  that  the  curve  Do  =  Do  (mo)  has  a  minimum  for  some  strictly 
positive  value  of  m 0.  For  large  mo  we  find  that  /i  ~  m q  and  hence  Do  ~  trq.  On  the 
other  hand,  numerical  computations  of  (A.l)  with  mo  =  0  yield  fi  =  0.1576654  >  0. 
It  follows  that  Do  ^  0.39707/mo  as  mo  — »  0+.  Thus  Do  blows  up  at  the  endpoints 
mo  — »■  0+  and  mo  00,  which  shows  that  this  curve  indeed  has  a  minimum. 
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